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Using a character expansion method, we calculate exactly the eigenvalue density of random 
matrices of the form M'M where M is a complex matrix drawn from a normalized distribution 
P(M) ~ exp(— Tr jAMBM'}) with A and B positive definite (square) matrices of arbitrary dimen- 
sions. Such so-called "correlated Wishart matrices" occur in many fields ranging from information 
theory to multivariate analysis. 



Physicists usually think of Wigner and Dyson as the 
fathers of random matrix theory [1,2]. However, twenty 
years before their first work on the subject, Wishart [3] 
examined random matrices of the form MM' as a tool 
for studying multivariate data. The properties of these 
so-called Wishart Matrices, which are viewed as "funda- 
mental to multivariate statistical analysis" [4], also find 
important applications in fields from information theory 
and communication [5-7] , to mesoscopics [8] , to high en- 
ergy physics [9], to econo-physics [10]. 

In many cases one is interested in Wishart matrices 
where the elements of M are not completely indepen- 
dent random variables, but have correlations along rows 
and/or columns. Important examples of this case occur 
in data analysis problems [11], where random samples 
have temporal and spatial correlations, and particularly 
in wireless communication and information theory [5,6]. 
The purpose of this paper is to derive the eigenvalue den- 
sity of correlated complex Wishart matrices exactly. A 
forthcoming longer paper will give more details of the 
derivation as well as discussing certain applications in 
depth. This problem has previously been studied in the 
limit of large matrices where perturbative expansions in 
1/A can be quite effective [11,6]. If either A or B is pro- 
portional to unity, simpler techniques can be used [12]. 

We first define the problem more precisely. Let M be 
an A by A' complex matrix chosen from a normalized 
distribution 

P(M) = 7r~ JVAr Vexp(-Tr {AMBM*}) (1) 

with A and B positive definite square matrices which 
define the correlations, and Tr indicates trace. Here, 
A/" -1 = detLA]^ detFB]^ and the factors of 7T, are nor- 
malization constants. An equivalent definition would be 
to let M = A _1 ' 2 ZB -1 ^ 2 where Z is a random complex 
matrix with independent entries of zero mean and unit 
covariance. Note that A is A by A and B is N' by N'. 
Without loss of generality, we can assume N > N' . For 
any operator O(M) we define the expectation bracket 
(O) to be an average over realizations of M so that 
(O) = JdM ■ 0(M)P(M). Note that the normalization 
is such that (1) = 1. 

Let A„ be the A' eigenvalues of M'M or equivalently 
the A' nonzero eigenvalues of MM' (we will also have 



A — A' eigenvalues of MM' precisely zero). Wc define 
the following quantities to calculate: 
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where 9 is the step function, A is assumed real, and in go- 
ing from Eq. 4 to 5 we have used lim c ^o I m log(— y+ie) = 
7r#(y) which is true for real y. The quantity we are most 
interested in is the eigenvalue density p(A). From Eqs. 
2-6 it is clear that we can obtain p by calculating G v (z). 
Below, we will show 

G v (z) = Q v (z)R v detLjj (7) 

Q^z)- 1 = A Ar ( a )A JV ,(fe)(-z) Ar '( JV '- 1 '/ 2 J l , (8) 

^=n£i>+*? (9) 

lu- llf ; V ':-V • , .,r v v ' (10) 

where i?„ is defined to be unity for A' > N — 1. In 
Eq. 7, Lij is an N by N matrix with elements Lij — 
g{a i b j - 1 v + A, z) for j < N' and L. tj = a\~ x for j > A' 
where we have defined a» and bj to be the eigenvalues of 
the matrices A and B. The function g is given by 

g(x;a,z) = x N - a e- zx T(a,-zx) (11) 

= // °° dA(A - z)"-^-** (12) 

with r the incomplete gamma function [14], and we note 
that for integer a > we have the simple form 

g(x, a, z) = x N - a (a - 1)! Y^a^T M (13) 
In Eq. 8 and throughout this paper we use the notation 
Av(*)= II (^•-^)=det[4" 1 ] (14) 

l<i<j<V 

to represent a T/-dimensional Vandermonde determinant. 
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From G u (z) we calculate C(A) using Eqs. 3 and 4. The 
differentiation (Eq. 3) with respect to v brings down a 
factor of log(A — z) in the argument of Eq. 12. This log 
becomes a step function when the limit is taken in Eq. 
4. We obtain 



JV 



C(X) = N' — Q (X)Ro det K , 
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where we also used G v —o (z) = 1. Here, we have defined 
TV by TV matrices with K<>f = gfabjiN, A) for n ± i 

for j > TV' and n ^ i. For 



and j < A' and KV- 
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the case of i = n we have Jf^ = e - a ^ x (N - 1)! for 

j < TV' and K { n f = for j > TV'. 

We then differentiate Eq. 15 (see Eq. 6) to obtain 
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where and T( nm ) are TV by TV matrices with ele- 

ments defined as follows: — for i ^ n and 



A' 



(n) 



[a n bj + (TV'(TV' - l)/2)/A]A'2 ) for j < TV' and 



K$> = for 7 > TV'. Also T ( " m) = A, ( ™ ) for i ^ m with 

T^ m) = (TV — l)a ro ^ g(a m bj,N - 1, A) for j < A' and 

7i7 } = for j > A'. 

The above expressions are our main results. The re- 
mainder of this paper comprises the proof of Eqs. 7-11 
from which all of our other results follow. We start by 
focusing on the case of square matrices M (so A = A'). 
We write 



Gu (z) = AAr-" 2 / dM e - Tr { AMBMt > det(MM+ 



we then define M = UmV where m is a diagonal ma- 
trix of the singular values m, of M and U and V are 
unitary matrices. We then separate the integral over M 
into integrals over the eigenvalues Aj = \rrii\ 2 of MM 1 
and "angular" integrals over U and V. This approach, 
common in random matrix theory [1], yields 

JV 



where U and V are A by TV unitary matrices which 
are integrated with the usual Haar measure over U (A) . 
Note that we have written the integral Da,b as a func- 
tion of A = mm^ (we will see this is indeed true). 
Here CA/"Aat(A) 2 Da,b(A) is precisely the joint probabil- 
ity density of the A's. As such, it is clear that this density 
(and also D) must vanish exponentially if any of the A's 
is taken to infinity (this will be important below). In 
particular, when A, B are identity matrices we can see 
that D a ,b(A) ~ exp(- J^i 

To address these integrals over U(N), we use the char- 
acter expansion method discussed in depth in Ref. [13]. 
This allows us to write 

e -T r {AU m VBVt ra tut} =ErarXr(AUmVBV , mtut) 

where a r arc expansion coefficients (discussed below) , the 
sum is over representations r of Gl(N), and \r is the 
character of the group element in the proper representa- 
tion. (Note that the representation theory of Gl(N) is 
identical to that of U(N)). A character is just the trace 
taken in the proper representation, so we have 

Xr (AUmVBvtmtut)=A r a6 c/>j/; e B^;; m r 9 c/ n " 

with lower repeated indices summed (and superscripts r 
tell us that the matrix is in representation r). We now 
use the orthogonality property [13] 
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with d r the dimension of representation r (discussed be- 
low). Combining the above three equations we obtain 



£>A, B (A) = Er^rd; 2 Xr(A)Xr(B)Xr(A) 



(18) 



As discussed in Ref. [13], each representation r is specified 
by a set of increasing integers < fcjv < fcjv-i • • • < fci, 
so the sum written over r is actually a ordered sum over 
the fc's (kj = A + rij — j in the notation of Ref. [13]). In 
Ref. [13] it it is also found that a r = s(k) det[l/(fcj + i — 

TV)!] = s(k)A N (k)/C(k) where C(k) = U?=iW~ Here 
s(k) = (-1)" withi; = A(A-l)/2-£\ kj. In the same 
work [13] it is also shown that a r /d r = s(fc)i 7 jv/C(fc) 



3 = 1 

where A = mm^ is the diagonal matrix of eigenvalues 
A, and J dX = YliLi Io° d\i and C is an A-dependent 
numerical constant (which we will not keep track of ex- 
plicitly but will fix at the end of the calculation). Here, 
Ajv(A) 2 is the Vandermonde determinant squared of the 
\ 's (which is the Jacobian of the transformation) and 



a r d~ 2 = s(k)F^/A N (k)C(k) 
The Weyl character formula tells us that [13] 
Xr (X) = det[x^}/A N (x) 



(19) 
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with kj the integers describing the representation r, and 
Xi the eigenvalues of X. Plugging Eq. 20 into Eq. 18 we 
obtain 
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A N (X) 2 D A , B (X) = A N (A)£ r $ r dct[A^] (21) 

$ r - a r d~ 2 Xr(A)xAB) (22) 
We next need a useful identity: 
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To show this we note that since A at (A) = Ili<j(^i — Xj) 
we can freely add a constant to each Aj and not change 



z). Next we use 



-I/saOIIJLi*?- 1 



Ajy. Thus, we have Ajv(A) = A^(A 
Ajv(a;i, . . . ,x N ) = A N (-l/xi, . . . , - 
so that we can relate A a? (A) to An(1/[z — A]). We then 
use l/(z — A) — 1/z = X/[z(z — X)] and we again shift each 
term in the Vandermonde determinant by —l/z. Finally 
we separate out factors of — z and write the Vandermonde 
determinant as on the far right of Eq. 14. 

Using the above expression for Aat(A) and plugging 
Eq. 21 into Eq. 17 yields 
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In Eq. 24 we have rewritten the determinants as sums 
over all permutations by using the completely antisym- 



metric Levi-Cevita tensor e r 



which is 1 if C\ , . 



,Cn 



is an even permutation of [1, . . . , N), is —1 if it is an odd 
permutation, and is otherwise zero. As mentioned above, 
the quantity (A^D) in the square brackets in Eq. 17 is 
exponentially convergent to zero when any Xi becomes 
large (Thus the quantity in the curly brackets is also ex- 



r di...d N i=l 



i=l ) 



(23) 
(24) 



ponentially convergent). Further, so long as a — 1 i= 
(which implies Cj — 1 + k^ ^ 0) the quantity in curly 
brackets goes to zero at the lower boundary Aj = 0. This 
enables us to trivially integrate by parts with respect to 
Aj where we differentiate the quantity in the curly brack- 
ets and integrate the quantity outside of the curly brack- 
ets and we do not obtain any boundary terms. We choose 
to do this integration exactly c, — 1 times to obtain 



G »w = (_,w(Li) /2 E ^.., N dAn(A«-zr +Ar - i (-ir- i E^ E e "lK di n 

Ci,...,cn " i=l r di...d N i=l p=l 



We would now like to interchange the order of inte- 
gration and summation such that all integrals are done 
first. However, if we did this we would end up with di- 
vergent integrals. To fix this problem, we insert a cutoff 



function such as /(A) = exp[— 5X] and at the end of the 
calculation we will take S to zero. (The precise form of 
the cutoff function will not matter). This allows us to 
reorder and write 



G "(*) - r _ Z \N{N -i)/2 H £c i- c « H fl(- 

■ ' r c\...c N d\...d N i=l 

We can now do the sums over c's and <f s to obtain 
G v {z) = N\ CM z- N{N - 1)/2 £ r $ r det[f#>] (25) 

(26) 



v + N -p 

P ^ = j™d\f(\)(\-zr+ N -i\ k > 
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The rather special form of the matrix expressed in Eq. 26 
allows us to calculate the determinant straightforwardly 
yielding det[i^ r) ] = A N (k) J" 1 . Thus we have 



$ r det[^ } ] 



ndet^jdet^int^P/r 
{-l) N ( N -V/ 2 A N {a)A N (b)J v 



where we have used Eq. 9, 19, 20, and 22 and the defini- 
tions of C(k) and s(k). Plugging this result into Eq. 



/ 4/(A,)(A,-,r-tn^ 



25 we now need only do the sum over r. This sum, 
as explained above, is actually a sum over < k^ < 
fc^v-i • • • <k\. Thus we have 



G v (z) = S E dct [a- 3 ]de% fe m1i^) (28) 

0<fejv...<fei 



(-1) 



k\ 



(29) 



with S = MCQ v {z) and we have absorbed the numerical 
constants Fn into C. We can now address the sum in Eq. 
28 using the Cauchy-Binet Theorem (see Appendix) to 
obtain G v (z) = S det[W (aibj)] with the function W(x) 
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defined by 

oo 

W(x) = ]T x k w{k) = J Q °° dX /(A) (A - z) 



v+N-l e -x\ 



k=0 



We can now remove the convergence function / (let- 
ting 5 — ► as discussed above) to obtain (See Eq. 12) 
W(x) — > g(x, v + N, z)x~ N . The factors of x~ N precisely 
cancel the prefactor AT and we recover the desired result 
Eq. 7 for N = N' (where R v = 1) up to the TV-dependent 
normalization prefactor C which we have not kept track 
of. To show that the normalization of Eq. 7 (i.e., C = 1) 
is indeed correct we need only verify Gq(z) = 1. To do 
this (using Eq. 13) we need to establish 



dct 
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m! 



Q( z ,oy 



which is easy to show using the Cauchy-Binet theorem 
(In Eq. 31, use w(k) = (-z) k {N - l)\/k\ for k < N - 1 
and w(k) — otherwise, so the determinants on the left 
hand side of Eq. 31 are precisely A (a) and A (6)) which 
completes the proof for the case of N = N'. 

Using the results we have derived for square matrices 
we can now easily derive results for rectangular matri- 
ces (N > N'). Given an N' dimensional matrix B with 
eigenvalues b\, . . . , b^' we consider an auxiliary N dimen- 
sional matrix B with the N' eigenvalues bi, . . . , 6jv' plus 
N — N' eigenvalues 6jv+i, ■ • ■ , &at- We then take a limit 
where 6jv'+i, bjy all go to infinity. By viewing the ma- 
trix M as being A _1 / 2 Z B -1 / 2 it is clear that taking this 
limit drives N — N' columns of M to zero and we obtain 
effectively an N by N' dimensional problem (with N — N' 
additional zero eigenvalues). To take these limits we will 
use the expansion (Ref. [14] Eq. 8.357) 



lim g(x, N + v, z) 



„N- 



N+v-\ 



N + v- 
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1 {N + v-l)(N + v-2) 



(30) 



As b?q — > oo we use the first term of this expansion and 
replace the j — N row of the matrix in Eq. 7 with 
(ai&Ar)^ 1 (— z) N+v ~ l . In the denominator (Q„) we have 
Aat(6) — > (&Ar) Ar_1 Aat_i(6) so the factors of b^ cancel 
to give a finite ratio. We next let b^-i — * oo. In taking 
this limit the first term in the expansion Eq. 30 would 
result in the j = N — 1 row of Lij being exactly pro- 
portional to the j = N row and thus we would obtain 
dct = 0. Thus, the leading divergence as 6jv-i — ► oo 
is actually from the second term of the expansion 30. 

row of L^ with 



N - 1 



We can then replace the j 

(N + v- l){a l b N - 1 ) N - 2 {-z) N+v - 2 . Again the diverg- 
ing powers of fe/v-i here are cancelled by powers in the 
denominator {Q v ) since A,/v-2(fr) = (&A r -i) A '~ 2 Ajv_i(6). 
This procedure can be continued until we have let all 
bx'+i ■ ■ .bjy — ► oo. We cancel all of the diverging 
terms then factor out the numerical prefactors (such as 
N + v — 1) to give R v , and factor out common factors 
of — z to obtain the general result quoted above in Eq. 
7-11 times (— z y( N ~ N ) which is due to the fact that, as 
mentioned above, our auxiliary problem has N — N' zero 
eigenvalues. In this way we complete our more general 
proof. 

Appendix: Cauchy-Binet Theorem [13]. Given TV- 
dimensional vectors a, and bi, and a function W(z) — 
Y^v=o w ^) zl convergent for \z\ < p then if \dibj\ < p for 
all i,j we have (where the determinants are all taken with 
respect to the indices i and j): 

J2 det[a^]det[b^]UtMki) = det[W(aibj)] (31) 

0<fc N ...<fci 
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